setwd("/Users/justingeorge/Library/CloudStorage/OneDrive-TexasA&MUniversity/Research/Defence Economics/JES_2025/Replication-materials")
df<-read.csv("tables_data.csv")
library(tidyverse)
library(ppcor)
library(purrr)

time<-c(1961,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015,2016,2017,2018,2019,2020,2021,2022,2023)

df<- df %>%
  filter(year %in% time)



spear_tables<-list()

# Compute Spearman correlations and p-values for each year
spear_tables[[1]] <- df %>%
  group_by(year) %>%
  summarise(
    spearman_corr = cor(gdp_2022, share, method = "spearman", use = "complete.obs"),
    p_value = cor.test(gdp_2022, share, method = "spearman", exact = FALSE)$p.value
  ) %>%
  mutate(
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    # Combine correlation, significance, and p-value into result_with_stars
    rho_12 = paste0(
      round(spearman_corr, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year, rho_12)  # Ensure dplyr::select is used

# Compute partial correlation by year
spear_tables[[2]] <- df %>%
  group_by(year) %>%
  summarise(
    # Compute partial correlation for each year
    partial_corr_result = list(
      pcor.test(
        x = gdp_2022,
        y = share,
        z = gdp_per_pop,
        method = "spearman"
      )
    ),
    .groups = "drop"  # Ungroup after summarisation
  ) %>%
  mutate(
    partial_corr_value = map_dbl(partial_corr_result, ~ .x$estimate),  # Extract correlation estimate
    p_value = map_dbl(partial_corr_result, ~ .x$p.value),             # Extract p-value
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    # Combine correlation, significance, and p-value into result_with_stars
    rho_123 = paste0(
      round(partial_corr_value, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year, rho_123)  # Ensure dplyr::select is used

# Compute partial correlation by year
spear_tables[[3]] <- df %>%
  group_by(year) %>%
  summarise(
    # Compute partial correlation for each year
    partial_corr_result = list(
      pcor.test(
        x = gdp_2022,
        y = share,
        z = c(gdp_per_pop, exp_border),
        method = "spearman"
      )
    ),
    .groups = "drop"  # Ungroup after summarisation
  ) %>%
  mutate(
    partial_corr_value = map_dbl(partial_corr_result, ~ .x$estimate),  # Extract correlation estimate
    p_value = map_dbl(partial_corr_result, ~ .x$p.value),             # Extract p-value
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    # Combine correlation, significance, and p-value into result_with_stars
    rho_1234 = paste0(
      round(partial_corr_value, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year, rho_1234)  # Ensure dplyr::select is used

# Save results to CSV
write.csv(spear_tables, "table_1.csv", row.names = FALSE)



df<-read.csv("tables_data.csv")

library(dplyr)

# Summarize the dataset
df_agg <- df %>%
  group_by(year) %>%
  summarise(
    gdp_total = sum(gdp_2022, na.rm = TRUE),
    pop_total = sum(pop, na.rm = TRUE),
    border_total = sum(exp_border, na.rm = TRUE),
    locattacks_total = sum(incidents_location, na.rm = TRUE),
    vicattacks_total = sum(incidents_target, na.rm = TRUE),
    me_total = sum(me,na.rm = TRUE),
    .groups = "drop"  # Ungroup the result after summarization
  )

df<- df  %>%
  left_join(df_agg, by="year")

df<- df %>%
  mutate(
    me_share = me/me_total,
    gdp_share = gdp_2022/gdp_total,
    pop_share = pop/pop_total,
    border_share = exp_border/border_total,
    ter_share_loc = incidents_location/locattacks_total,
    ter_share_vic = incidents_target/vicattacks_total
  )

df_wilcox <- df %>%
  mutate(
    avg1 = rowMeans(dplyr::select(., gdp_share, pop_share, border_share), na.rm = TRUE),
    avg2 = rowMeans(dplyr::select(., gdp_share, pop_share, border_share, ter_share_loc), na.rm = TRUE),
    avg3 = rowMeans(dplyr::select(., gdp_share, pop_share, border_share, ter_share_vic), na.rm = TRUE)
  )

time<-c(1990,1995,2000,2005,2010,2015,2020,2021,2022,2023)

df_wilcox<- df_wilcox %>%
  filter(year %in% time)

wilcox_table<-list()

wilcox_table[[1]] <- df_wilcox %>%
  group_by(year) %>%
  summarise(
    wilcox_stat = wilcox.test(me_share, avg1, paired = TRUE)$statistic,
    p_value = wilcox.test(me_share, avg1, paired = TRUE)$p.value
  ) %>%
  mutate(
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    share_1 = paste0(
      round(wilcox_stat, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year,  share_1)  # Ensure dplyr::select is used


wilcox_table[[2]] <- df_wilcox %>%
  group_by(year) %>%
  summarise(
    wilcox_stat = wilcox.test(me_share, avg2, paired = TRUE)$statistic,
    p_value = wilcox.test(me_share, avg2, paired = TRUE)$p.value
  ) %>%
  mutate(
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    share_2 = paste0(
      round(wilcox_stat, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year,  share_2)  # Ensure dplyr::select is used


wilcox_table[[3]] <- df_wilcox %>%
  group_by(year) %>%
  summarise(
    wilcox_stat = wilcox.test(me_share, avg3, paired = TRUE)$statistic,
    p_value = wilcox.test(me_share, avg3, paired = TRUE)$p.value
  ) %>%
  mutate(
    significance = case_when(
      p_value < 0.01 ~ "***",
      p_value < 0.05 ~ "**",
      p_value < 0.1  ~ "*",
      TRUE           ~ ""
    ),
    share_3 = paste0(
      round(wilcox_stat, 3),  # Correlation estimate
      significance,                 # Significance stars
      " (", round(p_value, 3), ")"  # P-value in parentheses
    )
  ) %>%
  dplyr::select(year,  share_3)  # Ensure dplyr::select is used

write.csv(wilcox_table, "table_2.csv", row.names = FALSE)
